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Abstract 

Developing  robust,  quantitative  methods  to  optimize  resource  allocations  in  response  to 
epidemics  has  the  potential  to  save  lives  and  minimize  health  care  costs.  In  this  paper,  we 
develop  and  apply  a  computationally  efficient  algorithm  that  enables  us  to  calculate  the 
complete  probability  distribution  for  the  final  epidemic  size  in  a  stochastic  Susceptible- 
Infected-Recovered  (SIR)  model.  Based  on  these  results,  we  determine  the  optimal  alloca¬ 
tions  of  a  limited  quantity  of  vaccine  between  two  non-interacting  populations.  We  compare 
the  stochastic  solution  to  results  obtained  for  the  traditional,  deterministic  SIR  model.  For  in¬ 
termediate  quantities  of  vaccine,  the  deterministic  model  is  a  poor  estimate  of  the  optimal 
strategy  for  the  more  realistic,  stochastic  case. 


Introduction 

As  rapid,  long-range  transportation  becomes  increasingly  accessible,  transmission  of  infectious 
diseases  is  a  growing  global  concern.  Advances  in  biomedical  therapies  and  production  have 
enabled  the  development  of  large  quantities  of  pre-pandemic  vaccine  [1],  The  United  King¬ 
dom,  Japan,  and  the  United  States  have  plans  to  stockpile  3.3  million,  10  million,  and 
40  million  doses,  respectively,  of  pre-pandemic  H5N1  vaccine  [2].  However,  in  the  face  of  a 
spreading  pandemic,  even  seemingly  extensive  resources  would  be  insufficient  to  provide  glob¬ 
al  coverage,  mandating  the  development  of  effective  protocols  for  the  allocation  of  limited 
vaccine  [3]  [4], 

A  starting  point  for  many  studies  of  disease  transmission  in  populations  is  the  Susceptible- 
Infected-Recovered  (SIR)  model  introduced  by  Kermack  and  McKendrick  [5],  In  this  model,  at 
any  given  time  each  individual  is  in  one  of  three  states.  The  dynamic  evolution  of  the  popula¬ 
tion  is  described  by  two  irreversible  transition  probabilities:  one  describes  the  rate  at  which  a 
susceptible  individual  becomes  infected,  and  the  other  describes  the  rate  at  which  an  infected 
individual  recovers  (or  dies).  The  net  effect  of  both  transition  rates  can  be  described  by  a  single 
number,  r0,  which  characterizes  how  effectively  the  infective  agent  moves  through  the  popula¬ 
tion.  On  average,  r0  describes  the  number  of  susceptible  individuals  infected  by  a  single  infected 
individual  in  a  population  of  susceptible  individuals.  Reducing  the  number  of  susceptible 
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individuals  in  a  population,  via  vaccination  for  example,  decreases  the  effective  reproductive 
number  re((.  When  reff  <  1,  the  number  of  infected  individuals  tends  to  decline,  whereas 
if  reff  >  1,  the  number  tends  to  grow. 

The  concept  of  a  reproductive  number  for  an  infectious  disease  can  be  generalized  to  more 
complex  models  of  epidemiology.  Previous  work  on  developing  optimal  vaccination  strategies 
typically  focus  on  minimizing  re(1,  either  by  proactive  dispersal  of  vaccine  before  the  infection 
reaches  a  population  [6],  or  reactive  dispersal  [7]  [8],  after  infection  has  been  detected  in  a 
group.  In  both  cases,  the  overall  size  of  the  epidemic,  measured  by  the  total  number  of  individ¬ 
uals  who  have  been  infected  throughout  the  course  of  the  epidemic,  is  lowered.  Vaccination  ac¬ 
complishes  this  by  removing  an  initial  number  of  susceptible  individuals  and,  thereby,  also 
suppressing  the  rate  of  infection. 

Numerous  computational  studies  of  large-scale  veterinary  infections,  such  as  foot-and- 
mouth  disease  [9]  [10],  Johne’s  Disease  [11],  as  well  as  human  infections  like  measles  [12]  and 
SARS  [13],  have  been  performed.  In  models  that  aim  to  capture  field  observations,  detailed, 
case-specific  information  such  as  demographics  of  the  population,  timing  and  logistics  for  vac¬ 
cine  deployment,  delays  associated  with  the  immune  response,  and  overall  vaccine  efficacy  are 
often  essential  to  the  investigation.  In  all  cases,  there  are  tradeoffs  between  complexity  and  real¬ 
ism,  and  between  computational  viability  and  the  generality  of  the  results. 

In  this  paper,  we  abstract  the  geographic,  demographic,  and  disease  specific  information, 
and  instead  focus  on  the  fundamental  problem  of  stochastic  SIR  dynamics  with  prophylactic 
vaccination  in  two  non-interacting  populations  (e.g.  two  well-separated  cities).  Previously, 
Keeling  and  Shattock  [2]  considered  the  deterministic  SIR  model  in  this  scenario,  and  obtained 
striking  results.  As  the  total  amount  of  available  vaccine  is  increased,  the  allocation  of  vaccine 
that  minimizes  the  total  number  of  infected  individuals  can  undergo  discontinuous  transitions. 
With  a  small  amount  of  vaccine,  the  optimal  strategy  involves  ensuring  that  the  smaller  popu¬ 
lation  was  well  protected  foremost.  However,  with  enough  vaccine,  the  optimal  strategy 
switches  abruptly  to  protecting  the  large  population,  leaving  the  smaller  population  entirely 
unprotected.  These  results  were  well  explained  in  terms  of  a  phenomenon  referred  to  as  “herd 
immunity,”  whereby  immunization  of  a  fraction  of  a  population  protects  even  those  who  are 
not  vaccinated,  by  reducing  the  effective  reproductive  number  reg  to  a  value  below  unity.  Vac¬ 
cination  removes  susceptible  individuals  from  the  population.  If  there  are  less  susceptible  indi¬ 
viduals  in  the  population,  on  average,  an  infected  individual  will  infect  fewer  individuals.  Herd 
immunity  occurs  when  rclT  <  1,  i.e.,  on  average,  at  the  start  of  the  epidemic  each  infected  indi¬ 
vidual  transmits  his  disease  to  less  than  one  person.  Keeling  and  Shattock  explained  the  sharp 
transitions  in  the  optimal  strategy  as  arising  from  a  strategy  that  aims  to  induce  herd  immunity 
in  the  largest  population  possible. 

While  the  deterministic  SIR  model  is  characterized  by  two  coupled  ordinary  differential 
equations,  the  stochastic  SIR  model  involves  a  high  dimensional  state  space  with  probabilistic 
transitions  between  partitions  of  the  overall  population,  characterized  by  the  number  of  indi¬ 
viduals  in  each  state.  Stochasticity  leads  to  noteworthy  differences  in  the  epidemic  size.  When 
re ff  >  1  the  probability  distribution  for  the  total  epidemic  size  is  bimodal  [14],  comprised  of  a 
roughly  Gaussian  peak  centered  at  the  deterministic  epidemic  size,  as  well  as  a  second  peak  for 
small,  “terminal  infections,”  describing  the  likelihood  the  disease  will  fail  to  propagate  signifi¬ 
cantly  during  the  initial  phase  of  infection.  The  peaks  of  the  distribution  are  well  separated 
when  the  population  size  is  large,  so  that  if  the  number  of  infected  individuals  exceeds  a  critical 
size,  the  epidemic  progresses  to  a  large  size,  characterized  on  average  by  the  deterministic  re¬ 
sults.  However,  the  non-negligible  probability  that  the  disease  will  fail  to  propagate  in  a  given 
population  results  in  significant  differences  for  optimal  allocation  of  vaccine  over  a  wide  range 
of  parameters. 
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The  rest  of  this  paper  is  organized  as  follows.  In  Methods  we  review  the  deterministic  SIR 
model,  and  its  stochastic  generalization.  We  approach  the  stochastic  problem  using  a  master 
equation  for  the  time  evolution  of  the  complete  probability  distribution  for  the  number  of  indi¬ 
viduals  in  each  state.  Building  on  the  computationally  efficient  algorithm  recently  developed  by 
Jenkinson  and  Goustias  [15],  we  introduce  a  modification  which  leads  to  even  greater  compu¬ 
tational  efficiency.  In  Results  we  compute  the  probability  distribution  for  the  final  epidemic 
size  for  a  range  of  parameters  to  identify  regimes  for  which  the  stochastic  and  deterministic 
models  differ  most  significantly.  We  compute  the  optimal  allocation  of  vaccine  between  two 
non-interacting  populations,  and  compare  our  results  with  the  deterministic  case.  Stochastic 
effects  are  most  pronounced  in  situations  involving  an  intermediate  amount  of  resource  avail¬ 
ability.  We  conclude  with  a  discussion  of  our  results  and  future  directions. 


Methods 

We  briefly  review  the  deterministic  SIR  model  [5],  a  system  of  coupled  differential  equations 
for  modeling  the  growth  of  an  epidemic  in  a  well-mixed  population  within  which  all  agents  in¬ 
teract  equally  with  all  other  agents.  The  model  describes  a  population  of  N  individuals  divided 
into  three  classes,  susceptible  S,  infected  I,  and  recovered  R: 

^  =  -/fS(*)i(0;  (i) 

®  =  /?s(t)i(0-yi(0;  (2) 


dR{t) 

dt 


(3) 


We  may  omit  the  equation  for  the  recovered  class  because  we  can  always  deduce  the  number  of 
recovered  individuals  from  the  fact  that  the  total  number  of  individuals  in  the  population,  N,  is 
fixed,  so  that  R(f)  =  N-S(t)-I(t). 

Equations  1-3  can  be  thought  of  as  a  mean  field  theory  where  the  continuous  variables  S(f) 
and  1(f)  are  the  average  values  (over  many  iterations)  of  two  discrete  integer-valued  variables  S 
and  I.  At  any  time  then,  the  system  can  be  characterized  as  being  in  a  state  (S,  f),  which  can  un¬ 
dergo  one  of  two  transitions: 

(S,  7)  -►  (S-  1,1+  l)at  rat efSSI; 

(S,  I)  — >•  (S,  I  —  1)  at  rate  yl. 


The  parameters  / 3  and  y  can  be  defined  in  terms  of  physical  observables,  the  average  number  of 
contacts  each  person  makes  per  day  c,  the  probability  of  infection  through  contact  p,  and  the 
characteristic  duration  of  the  infection  T: 


P  = 


f  rate  at  which  each  \  f  probability  of  infection  \ 
y  individual  makes  contacts  J  y  from  contact  J  exp 
( total  population  size  N)  N 


(4) 


y  = 


i 

characteristic  duration  of  infection 


1 

T' 


(5) 
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For  each  set  of  parameters  /i  and  y,  the  reproductive  number,  r0  is  defined  to  be: 


B  S„  S„ 

ro=  — =  cxpxTx-, 


(6) 


where  S0  is  the  initial  number  of  susceptible  individuals,  S(f  =  0).  Here  r0  can  be  interpreted  as 
the  average  number  of  new  infections  a  single  infected  individual  will  produce  in  a  completely 
susceptible  population.  Thus  if  r0  <  1,  in  the  deterministic  model  dl/dt  <  0,  the  number  of  in¬ 
fected  individuals  will  decline  from  the  initial  seed  value,  1(f)  <  I0  =  I(f  =  0),  and  no  epidemic 
will  occur.  In  our  numerical  simulations,  the  value  of  r0  is  tuned  by  varying  fi  Because  ft  is  in¬ 
versely  proportional  to  N,  transmission  is  frequency  dependent,  and  the  rate  at  which  each  in¬ 
dividual  makes  contacts  with  others,  c,  is  independent  of  the  population  size  N. 

In  this  paper  we  investigate  the  effects  of  prophylactic  vaccination.  Vaccinating  V  individu¬ 
als  proactively  corresponds  to  removing  V  susceptible  individuals  before  the  epidemic  begins, 
thus  lowering  the  effective  reproductive  number.  This  assumes  that  the  vaccine  is  completely 
effective.  Let  r0  denote  the  reproductive  number  prior  to  vaccination,  and  re((  denote  the  effec¬ 
tive  reproductive  number  which  is  achieved  after  vaccinating  V  individuals: 

^  SH  -  V  /Sn  -  V\ 

,a=«Frx  — =  (7) 


If  a  sufficient  number  of  individuals  V  are  vaccinated,  rdf-  may  be  reduced  to  a  value  below 
unity,  so  that  dl/dt  <  0,  and,  as  a  result,  the  epidemic  will  not  grow.  Thus  the  entire  population 
will  be  safeguarded  without  vaccinating  the  entire  population.  This  phenomenon  is  known  as 
herd  immunity. 

In  the  stochastic  SIR  model,  the  infection  and  recovery  reactions  are  modeled  as  continu¬ 
ous-time  Markovian  processes.  Let  be  the  probability  at  time  t  of  a  population  with  S  sus¬ 
ceptible  individuals  and  I  infected  individuals  with  N  =  S0+I0  where  S0  and  I0  are  the  initial 
values  of  S  and  I.  The  evolution  of  in  time  is  then  governed  by: 

=  P(S+  X)(J  ~  !)  x  0s+m— i(0 

+  H-f  + 1)  x  ^s,/+ iM  (8) 

-  (psi  +  yl)  x  <j)SJ(t), 


where  the  first  two  terms  on  the  right  hand  side  correspond  to  transitions  into  the  state  (S,  I) 
by  a  susceptible  individual  becoming  infected  or  an  infected  individual  recovering,  respectively, 
and  the  third  term  corresponds  to  the  probability  of  leaving  the  state  (S,  I)  through  infection  or 
recovery.  While  the  deterministic  model  tracks  the  time  evolution  of  two  ensemble  averaged 
variables  S(f)  and  1(f),  the  stochastic  model  has  up  to  (S0  +  f0  +  l)(So  +  1)  ~  N2  possible  states 
(. I  can  take  values  0  to  S0  +  Iq  and  S  can  take  values  0  to  S0).  All  the  probabilities  <f>sj(t)  can  be 
assembled  into  a  vector  (f>  and  the  entire  system  of  equations  can  then  be  written  in 
matrix  form. 

We  computationally  integrate  this  system  of  equations  using  a  modified  version  of  Jenkin- 
son  and  Goutsias’s  method  [15]  of  Implicit-Euler  integration.  The  matrix  A  consists  of  the  co¬ 
efficients  of  Equation  8,  and  describes  the  transition  probabilities: 

dtj>  =  A  (f)dt .  (9) 
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The  above  equation  is  discretized  by  introducing  a  time  step,  which  controls  the  accuracy  of 
the  method: 


4>ti+1  =  4>u  +  A^i+1(f.+i  -  0; 

=  (i- aa  ty%,. 


(10) 


The  ordering  of  the  components  of  the  vector  (6  in  such  a  way  so  that  A  is  lower  triangular  re¬ 
duces  the  number  of  computations  needed  to  solve  Equation  10  from  0(KS)  to  0(K2),  where 
K  is  the  length  of  the  vector  0  and  scales  with  the  system  size.  [15]  We  made  modifications  to 
the  way  the  algorithm  counts  states,  enabling  considerably  faster  computational  speed,  espe¬ 
cially  as  the  population  size  increases.  Where  Jenkinson  and  Goutsias  take  the  approach  of 
counting  the  so-called  “degree  of  advancement,”  a  scenario  in  which  each  state  corresponds  to 
a  specific  sequence  of  reactions,  we  instead  take  the  “population  process”  approach  by  enumer¬ 
ating  all  states  of  the  system  without  tracking  which  reactions  might  have  led  the  system  to  the 
state  in  question.  In  both  methods  one  begins  with  {S0  + 10+  1)(S0  +  1)  states.  In  our  method, 
we  remove  those  states  that  have  zero  probability  of  occurring,  but  are  included  in  the  original 
computational  algorithm.  For  example,  many  states  where  S  +  I>  N  are  retained  in  the  degree 
of  advancement  procedure  but  are  explicitly  excluded  in  our  method.  The  result  is  that  we 
track  [(S0  +  l)(/0  +  1)  +  (S0  +  l)(S0)/2]  states,  which  in  the  limit  I0  <C  S0  is  approximately 
~  N2 12.  As  the  system  size  N  grows,  the  difference  in  the  total  number  of  states  between  the 
two  methods  can  significantly  impact  the  time  it  takes  to  integrate  the  system  of  equations. 

The  system  is  initialized  with  a  population  size  N,  I0  infected  individuals,  and  an  initial  re¬ 
productive  number  r0.  Thus,  at  time  t-  0  the  probability  of  state  (f>N  _  j  j  ( 0)  =  1  and  the  proba¬ 
bility  of  all  other  states  equals  zero.  The  collection  of  probabilities  6  of  all  accessible  states  is 
then  evolved  forward  in  time  until  the  distribution  reaches  a  stationary  state  where  the  proba¬ 
bility  of  having  any  state  (S,  I)  where  I  >  0  is  vanishingly  small.  At  that  point,  all  individuals  in 
the  initial  population  of  size  N,  have  either  been  infected,  and  are  now  recovered,  or  remain 
susceptible.  For  the  parameters  considered  here,  we  observe  that  an  integration  time  of  t  =  200 
is  sufficient  in  all  cases.  Once  the  simulation  is  complete,  we  define  the  final  epidemic  size  E  as: 


E  =  limJV  —  S(f)in  the  stochastic  model; 

^°°  (11) 

E  =  lim  N  —  S(f)  in  the  deterministic  model; 

t — »oo 


Fig.  1  illustrates  numerical  results  for  the  epidemic  size  distribution  P(E),  describing  the 
probability  of  having  a  total  of  E  individuals  infected  over  the  course  of  the  entire  simulation 
period.  We  observe  that  for  r0>  1  the  probability  distribution  consists  of  two  parts.  On  the  left 
side  of  Fig.  1  A,  there  is  peak  describing  small,  “terminal  infections,”  which  fail  to  propagate  sig¬ 
nificantly  in  the  population  (i.e.  the  infection  terminates  before  a  large  number  of  individuals 
are  impacted).  The  peak  describing  terminal  infections  decays  approximately  exponentially 
from  the  peak  value  at  P(/0).  In  the  stochastic  model,  there  is  always  a  nonzero  probability  the 
infection  will  end  without  becoming  a  large  scale  epidemic. 

When  r0  >  1,  the  distribution  P(E)  exhibits  a  second  peak  towards  the  right  side  of  Fig.  1A, 
describing  “large-scale  epidemics.”  This  peak  is  approximately  centered  at  the  epidemic  size 
predicted  by  the  deterministic  SIR  model,  illustrated  for  each  value  of  r0  by  the  corresponding 
vertical  dashed  line  in  Fig.  1  A.  The  size  of  the  large-scale  epidemic  scales  with  the  size  of  the 
population  N,  resulting  in  increasing  separation  of  the  peaks  for  increasing  population  sizes. 

To  quantify  the  likelihood  of  a  terminal  infection  versus  a  large-scale  epidemic  by  the  relative 
weight  associated  with  each  of  the  peaks,  numerically  we  define  the  point  separating  the 
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Fig  1.  Epidemic  size  distribution.  Figure  A  illustrates  the  final  epidemic  size  distributions  P(E)  for  various  values  of  r0  in  the  stochastic  SIR  model.  The 
corresponding  result  for  the  deterministic  SIR  model  (a  number  E)  is  marked  in  each  case  by  a  dashed  line.  The  left  vertical  scale  in  Figure  A  describes  the 
small  terminal  infections  (left  of  the  first  horizontal  scale  break),  and  the  right  vertical  scale  describes  large-scale  infections  (right  of  the  first  horizontal  scale 
break).  Figure  B  illustrates  the  corresponding  cumulative  probability  distribution  P(E  <  Emax)  or  the  probability  of  an  epidemic  of  size  less  than  Emax.  This  also 
shows  the  relative  weight  in  the  terminal  infection  and  in  the  large-scale  epidemic.  In  each  case  N  =  500  and  /0  =  1 . 

doi:10.1371/journal.pone.0115826.g001 


terminal  infection  and  the  large-scale  epidemic  to  correspond  to  the  local  minimum  in  proba¬ 
bility  that  exists  between  the  two  peaks.  The  likelihood  of  a  terminal  infection,  represented  by 
the  total  weight  in  the  terminal  infection  peak,  decreases  with  increasing  values  of  r0  and  70.  As 
r0  approaches  unity  from  above,  the  large-scale  epidemic  progressively  decreases  in  mean  size, 
but  increases  in  variance.  Eventually  the  distinction  between  terminal  infections  and  large- 
scale  epidemics  vanishes  (the  local  minimum  in  P(E)  ceases  to  exist).  This  is  associated  with  a 
critical  phase  transition  [16],  and  occurs  at  a  value  of  r0  that  approaches  unity  as  the  population 
size  N  tends  to  infinity.  When  r0<  1,  the  probability  distribution  is  described  only  by 
terminal  infections. 

The  cumulative  epidemic  size  distribution  P(E  <  fjmax)  describes  the  probability  of  having 
an  epidemic  of  size  less  than  £max,  and  is  shown  in  Fig.  IB.  The  extended  flat  portions  of  the 
curves  indicate  that  a  population  of  N  =  500  individuals  is  well  within  the  large  population 
limit,  defined  by  a  large  separation  between  the  terminal  infection  and  large-scale  epidemic, 
with  little  probability  of  observing  an  epidemic  size  in  between  the  two.  Fig.  IB  also  illustrates 
how  the  total  probability  is  distributed  between  the  terminal  infection  and  the  large-scale 
epidemic.  The  smaller  the  value  of  r0,  the  greater  the  likelihood  that  the  initial  seed  population 
of  infected  individuals  will  fail  to  spread  the  disease. 


Results 

Our  aim  is  to  highlight  key  differences  between  stochastic  and  deterministic  approaches  to  de¬ 
veloping  a  framework  for  the  optimal  allocation  of  vaccine  between  two  non-interacting  popu¬ 
lations.  We  begin  by  observing  how  a  single  population  reacts  to  different  levels  of  vaccination 
in  the  stochastic  and  deterministic  SIR  models;  the  results  provide  the  basis  for  the  optimiza¬ 
tion  process.  Subsequently,  we  determine  the  optimal  deterministic  and  stochastic  solutions 
that  minimize  the  average  epidemic  size  between  two  populations,  where  one  is  twice  the  size 
of  the  other,  and  contrast  their  properties.  We  then  demonstrate  the  robustness  of  these  results 
by  computing  the  corresponding  optimal  solutions  over  a  range  of  alternative  parameters,  that 
include  variations  in  the  ratio  of  population  sizes,  and  increasing  the  number  of  individuals 
who  are  initially  infected.  Finally,  we  consider  an  alternative  optimization  based  on  imposing  a 
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maximum  tolerance  for  the  epidemic  size  and  show  that  the  stochastic  optimal  solution  better 
fulfills  this  measure  than  the  deterministic  optimum. 


Impact  of  Vaccination  on  the  Epidemic  Size  of  a  Single  Population 

We  first  consider  how  the  epidemic  size  within  a  single  population  decreases  as  a  function  of 
increasing  vaccine  allocation.  Vaccine  allocation  V  removes  V  susceptible  individuals  from  the 
initial  state  (S0,  Iq)  — >  (So  -  V,  I0 )  after  which  the  stochastic  SIR  model  evolves  according  to 
Equation  8  (Equations  1-3  for  the  deterministic  SIR  model).  The  resulting  dynamics  determine 
the  size  of  the  epidemic  according  to  Equation  1 1 .  Decreasing  the  initial  number  of  susceptible 
individuals  S0  by  V  will  not  in  general  lead  to  a  corresponding  reduction  V  in  the  final  epidemic 


A 


C 


Number  Vaccinated  (V) 


Fig  2.  Comparison  of  the  effects  of  vaccination  between  stochastic  and  deterministic  models  for  various  values  of  r0and  l0.  Figures  A,  B,  and  C 

illustrate  the  epidemic  size  as  a  function  of  the  number  of  individuals  vaccinated  V.  For  each  value  of  r0,  the  solid  lines  shows  the  average  epidemic  size  (E)  = 
/  P(E)  x  E  dE,  and  dashed  lines  represent  the  corresponding  deterministic  epidemic  size  E.  The  deterministic  solution  does  not  depend  significantly  on  /<,  and 
is  computed  with  /0  =  1 .  Figure  D  illustrates  the  standard  deviation  in  the  large-scale  epidemic  for  various  values  of  r0.  Results  are  independent  of  l0  for  the 
range  of  l0  shown.  While  the  average  epidemic  size  decreases  monotonically  with  V  (Figures  A-C),  the  standard  deviation  of  the  large-scale  epidemic  size 
initially  increases  with  increasing  V  (decreasing  reff),  reflecting  greater  variability  as  reff  is  reduced,  similar  to  the  increase  in  the  width  of  the  large-scale 
epidemic  peak  with  decreasing  r0  in  Fig.  1  A.  The  results  obtained  here  are  for  a  population  of  N  =  500. 


doi :  1 0. 1 371  /jou  rnal.pone.01 1 5826  .g002 
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size.  An  important  quantity  for  optimizing  the  allocation  is  the  incremental  reduction  in  the 
expected  epidemic  size  per  incremental  increase  in  the  allocation. 

In  Fig.  2  we  illustrate  the  numerical  results  for  a  population  of  N  -  500  individuals  with  dif¬ 
ferent  initial  numbers  of  infected  individuals  I0  and  different  reproductive  numbers  r0.  An 
amount  of  vaccine  V  (0  <  V  <  N  - 10)  is  given  to  the  population  and  we  compute  the  average 
final  epidemic  size  (E)  =  f  P{E)  x  E  dE  as  a  function  of  V,  where  P(E)  is  computed  as  in 
Fig.  1  A.  We  also  plot  the  corresponding  deterministic  curve  in  each  case,  where  P(E)  here  is  de¬ 
scribed  by  a  d-function  at  the  deterministic  epidemic  size  (5(E).  In  the  stochastic  model,  the 
quantity  (E)  depends  on  the  statistics  of  both  the  terminal  infection  and  the  large-scale  epi¬ 
demic;  (E)  may  not  correspond  to  an  epidemic  size  that  is  likely  to  be  observed  because  there 
may  be  a  large  separation  between  the  observed  sizes  of  terminal  infections  and  large-scale  epi¬ 
demics,  with  the  mean  size  lying  somewhere  in  between. 

Flerd  immunity  occurs  in  the  deterministic  SIR  model  when  the  initial  effective  growth  rate 
of  the  number  of  infected  individuals  in  the  population  becomes  less  than  unity  (reff  <  1)  [17], 
and  is  achieved  at  a  value  of  V  determined  by  Equation  7,  i.e.  when  V7S0  =  V/(N  -  I0)  =  1  - 
r0_1.  In  the  limit  of  large  populations,  the  fraction  that  must  be  vaccinated  to  achieve  herd  im¬ 
munity  approaches  1  -  r0~l.  Thus  for  a  population  N  =  500,  herd  immunity  occurs  approxi¬ 
mately  when  V  =  250  for  r0  =  2,V=  400  for  r0  =  5,  and  V  =  450  for  r0=  10.  Approaching  this 
value,  the  incremental  reduction  in  expected  epidemic  size  per  increase  in  vaccine  allocation 
increases  monotonically.  Note  that  the  peak  epidemic  reduction  rate  occurs  for  a  slightly 
smaller  V  when  N  is  finite,  compared  to  the  N  — >  oo  limit.  This  is  due  to  the  non-negligible 
(compared  to  N)  contribution  of  the  initial  seed  population  of  infected  individuals  I0  in  the  def¬ 
inition  of  the  herd  immunity  threshold  V=  (1  -  r0-1)(N-  70). 

In  the  stochastic  model,  the  corresponding  transition  is  subtler.  Increasing  the  vaccine  allo¬ 
cation  has  three  effects  on  P(E).  It  decreases  the  mean  size  (E)  and  increases  the  variance  of  the 
large-scale  epidemic,  and  also  increases  the  relative  likelihood  of  terminal  infections.  We  asso¬ 
ciate  the  onset  of  “effective  herd  immunity”  in  the  stochastic  model  with  the  value  of  V  for 
which  the  distinction  between  terminal  infections  and  large-scale  epidemics  ceases  to  exist,  as 
measured  by  the  existence  of  a  local  minimum  in  P(E).  Because  of  the  probability  of  terminal 
infections,  this  generally  occurs  at  a  value  of  V  which  is  smaller  than  that  of  the  herd  immunity 
transition  in  the  deterministic  model.  Furthermore,  unlike  the  deterministic  model,  in  the  sto¬ 
chastic  case  approaching  the  onset  of  effective  herd  immunity  does  not  coincide  with  a  specific 
value  of  reff  and  is  not  generally  the  point  of  maximum  impact  per  vaccine  in  the  allocation  (as 
measured  by  reduction  in  the  average  epidemic  size). 

Before  the  population  has  reached  the  deterministic  herd  immunity  transition,  i.e.  when  refj 
>  1,  the  deterministic  epidemic  size  E  (dashed  lines  of  Fig.  2A,  B,  C)  is  generally  larger  than  the 
average  stochastic  epidemic  size  (E)  (solid  lines).  While  the  maximum  size  of  the  large-scale 
epidemic  can  be  greater  than  the  deterministic  epidemic  E,  the  average  epidemic  size  (E)  is 
smaller  due  to  the  fact  that  the  stochastic  model  includes  the  possibility  of  a  terminal  infection. 

When  sufficient  vaccine  is  available  to  establish  herd  immunity  in  the  deterministic  model, 
the  situation  is  reversed,  and  the  deterministic  size  E  is  generally  smaller  than  the  average  sto¬ 
chastic  epidemic  size  ( E ).  When  re^  <  1,  in  the  deterministic  model  dlldt  <  0  at  t  =  0,  and  the 
initial  number  of  infected  individuals  decreases.  On  the  other  hand,  stochastically  there  is  al¬ 
ways  a  possibility  that  the  initial  number  of  infected  individuals  will  grow.  Hence,  beyond  the 
herd  immunity  threshold,  the  deterministic  epidemic  E  is  smaller  than  the  average  stochastic 
epidemic  (E). 

The  size  of  the  average  stochastic  epidemic  ( E )  also  approaches  the  deterministic  outcome  E 
as  both  I0  and  r0  become  large.  A  larger  value  of  r0  causes  each  infected  individual  to  infect 
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more  susceptible  individuals,  while  larger  values  of  I0  makes  it  less  likely  for  every  member  of 
the  initial  group  of  infected  individuals  to  recover  before  spreading  their  disease.  Both  of  these 
effects  decrease  the  probability  of  a  terminal  infection. 

Fig.  2D  shows  the  standard  deviation  of  the  large-scale  epidemic.  While  this  quantity  does 
not  factor  independently  into  any  of  the  optimization  problems  considered  in  this  paper,  the 
variation  of  the  standard  deviation  of  the  large-scale  epidemic  with  r0, 10,  and  V  illustrates  sev¬ 
eral  key  features  of  the  stochastic  model  that  differ  from  the  deterministic  case  where  the  stan¬ 
dard  deviation  is  a  priori  zero. 

Firstly,  the  standard  deviation  for  large  epidemics  is  independent  of  the  initial  number  of  in¬ 
fected  individuals  (i.e.  I0  =  1,2,  5),  as  long  as  /<>  <C  N.  If  the  infection  grows  into  a  large-scale  ep¬ 
idemic,  reaching  a  size  comparable  to  the  system  size  N,  the  impact  of  the  original  number  of 
initial  infected  individuals  I0  on  the  SIR  dynamics  becomes  negligible.  The  standard  deviation 
does,  however,  increase  with  decreasing  r0.  This  is  illustrated  directly  for  several  values  of  r0  in 
Fig.  1  A.  In  each  case,  a  smaller  value  of  r0  implies  each  infected  individual  on  average  infects 
fewer  susceptible  individuals.  This  increases  the  variability  in  the  outcome  arising  from 
stochastic  effects. 

The  standard  deviation  in  the  large-scale  epidemic  is  also  a  function  of  V,  the  amount  of 
vaccine  allocated  to  the  population.  When  V  is  small,  the  standard  deviation  increases  with  in¬ 
creasing  V.  This  is  attributed  to  the  fact  that  having  more  vaccinated  individuals  reduces  re(T. 
This  effect  is  balanced  by  the  fact  that  more  vaccinated  individuals  results  in  fewer  available 
configurations  (S,  I)  for  the  system  to  transition  into.  Hence  as  V  increases  further,  the  stan¬ 
dard  deviation  eventually  peaks  and  then  drops  sharply  to  zero.  The  value  of  zero  corresponds 
to  the  disappearance  of  the  local  minimum  in  P(E)  separating  terminal  infections  from  large- 
scale  epidemics,  coinciding  with  our  definition  of  effective  herd  immunity.  The  observation 
that  the  standard  deviation  peaks  at  a  value  of  V  just  below  the  onset  of  effective  herd  immunity 
indicates  that  the  largest  uncertainty  in  the  size  of  the  large-scale  epidemic  is  expected  for  allo¬ 
cations  just  below  effective  herd  immunity. 

An  important  quantity  in  determining  the  optimal  allocation  of  vaccine  is  the  “gain”  G, 
which  corresponds  to  the  incremental  reduction  in  the  expected  epidemic  size  (E)  per  incre¬ 
mental  increase  in  the  allocation  V: 


i  /  ~u\ 

G  = - in  the  stochastic  model; 

dv 

dE 

G  =  — —  in  the  deterministic  model; 
dV 


(12) 


In  the  stochastic  model,  the  dependence  of  the  gain  on  r0  can  be  separated  into  three  distinct 
cases:  the  subcritical  case  r0  <  1  (not  shown),  the  large  r0  case  (r0  a  2.5  for  the  other  parame¬ 
ters  considered  here),  and  the  intermediate  case  where  r0  is  greater  than  the  critical  value  of 
unity,  but  below  the  large  r0  limit  (1  <  r0  H  2.5).  For  r0  S  2.5,  the  gain  initially  increases  (at  a 
smaller  rate  than  the  corresponding  deterministic  curve),  peaks,  and  then  declines  to  zero.  In 
this  case  then,  when  optimizing  vaccine  allocation,  there  is  a  value  of  V  prior  to  reaching  herd 
immunity  where  the  gain  from  vaccination  peaks.  This  is  shown  for  r0  =  5  in  Fig.  3A.  For  1  < 
r0  5;  2.5,  the  gain  G  instead  declines  continuously  from  a  maximum  value  at  V  =  1.  This  is  illus¬ 
trated  explicitly  for  r0  =  2  in  Fig.  3B.  This  behavior  implies  that  the  larger  the  vaccine  allocation 
V  given  to  the  population,  the  smaller  the  benefits  of  even  larger  V.  When  r0<  1  (not  shown), 
there  is  essentially  zero  probability  of  a  large-scale  epidemic,  thus  very  little  decrease  in  the  epi¬ 
demic  size  (E)  per  increase  in  allocation  V,  and  the  gain  is  effectively  zero. 

In  the  deterministic  model,  for  all  r0>  1  the  gain  curve  follows  the  same  qualitative  pattern 
as  the  large  r0  case  in  the  stochastic  model.  Initially,  the  gain  G  increases,  rising  sharply  prior  to 
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Number  Vaccinated  (V)  Number  Vaccinated  (V) 


Fig  3.  Gain  G  as  a  function  of  the  vaccine  allocation  V.  For  both  the  stochastic  (solid  lines)  and  deterministic  (dashed  lines)  models,  Figs.  3A  and  3B 
illustrate  G  (i.e.  slopes)  of  the/0  =  1  curves  in  Figs.  2A  and  2B,  respectively.  Figure  A  shows  G  for  r0  =  5  for  both  a  population  of  N  =  500  individuals  and  also 
one  of  N  =  1 000  individuals.  Figure  B  shows  the  corresponding  results  for  r0  =  2.  The  gain  curves  for  a  population  of  N  =  500  terminate  aXV=  500  and  those 
for  a  population  of  N  =  1 000  terminate  at  V=  1 000. 

doi:  1 0. 1 371  /journal  .pone.0 1 1 5826.g003 


herd  immunity,  and  then  falling  sharply  after  the  vaccine  exceeds  the  herd  immunity  point. 
This  means  there  is  a  significant  increase  in  the  gain  G  from  vaccination  as  the  level  of  vaccine 
in  the  population  nears  the  herd  immunity  threshold. 

What  is  notably  different  between  the  stochastic  and  deterministic  models  is  the  sharpness 
of  the  peak  and  the  rate  of  decline  that  follows.  This  is  apparent  when  comparing  the  stochastic 
and  deterministic  curves  of  Fig.  3  A.  This  illustrates  that  beyond  a  threshold  level  of  vaccination 
{reg  <  1  in  the  deterministic  model)  there  is  almost  no  reduction  in  E  by  increasing  the  alloca¬ 
tion  V  to  the  population.  In  the  stochastic  model  there  is  not  as  definitive  a  threshold  level  of 
vaccination.  The  gain  G  in  the  stochastic  model  begins  to  decline  well  before  effective  herd  im¬ 
munity  is  reached.  Thus  in  the  stochastic  model,  the  point  of  diminishing  returns  from  vacci¬ 
nation  will  generally  take  place  at  smaller  vaccine  allocations  V  compared  to  the 
deterministic  model. 


Optimal  Vaccination  Allocation  for  Two  Populations 

Next  we  consider  the  problem  of  vaccine  allocations  for  two  non-interacting  populations  (e.g., 
two  well-separated  cities).  This  scenario  isolates  a  fundamental  tradeoff  in  resource  manage¬ 
ment,  whereby  allocating  vaccine  to  one  population  occurs  at  the  expense  of  the  other. 

Unless  otherwise  specified,  we  identify  properties  specific  to  each  population  with  super¬ 
scripts  1  and  2.  We  assume  in  this  section  one  population  is  relatively  small  ( N 1  =  500  individu¬ 
als),  and  the  other  is  relatively  large  ( N 2  =  1000  individuals).  Both  populations  are  initialized 
with  a  single  infected  individual  I*  =  I2  =  1.  In  this  scenario,  both  are  also  exposed  to  an  infec¬ 
tion  with  the  same  reproductive  number  r'}  =  r%,  so  that  /31  =  2/?2,  as  per  Equation  4.  A  fixed 
total  amount  of  vaccine  V  (0  <  V  <  1498,  where  the  maximum  value  of  V  is  given  by  N1  + 

N2  —  1^  —  =  1498,  accounting  for  one  seed  infected  individual  in  each  population)  can  be 

distributed  between  the  two  populations,  so  that  the  small  population  receives  V1,  and  the  large 
population  receives  V2  =  V-  V1.  We  define  the  optimal  allocation  to  be  the  partition  of  V  into 
[V1,  V2]  =  [V1,  V-  U1]  that  minimizes  the  average  total  final  epidemic  size  (E)  =  (E1)  +  (E2), 
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Fig  4.  Optimal  Solution  in  the  Stochastic  Model.  The  black  lines  of  Figures  A  and  B  illustrate  the  optimal  allocation  of  an  amount  of  vaccine  t/that 
minimizes  (£}  for,  respectively,  r0  =  5  and  r0  =  2.  The  y-axis  here  shows  the  amount  of  vaccine  allocated  to  the  smaller  population  I/1 .  The  color  scale 
indicates  the  average  epidemic  size  (E)  corresponding  to  a  particular  allocation  of  vaccine.  Figure  C  shows  how  the  optimal  solution  varies  as  r0  changes 
from  2  to  5.  Switching  behavior  vanishes  when  r0  rs  2.9  (not  shown).  Results  are  obtained  for  two  populations  of  sizes  A/1  =  500  and  N2  =  1 000.  Both 
populations  are  initialized  with  a  single  infected  individual  l'0  =  P0  =  1 . 
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where  0  <  E<  N1  +  N2  =  N  =  1500.  In  this  scenario,  the  cost  of  producing  and  distributing  vac¬ 
cine  is  not  taken  into  account,  so  it  is  always  beneficial  to  use  all  of  the  available  vaccine. 

Our  objective  is  to  determine  the  optimal  solution  as  a  function  of  V for  the  stochastic  SIR 
model.  We  compare  our  results  to  the  corresponding  optimal  solution  for  the  deterministic 
SIR  model,  which  we  also  compute.  For  both  models,  we  sample  the  space  of  all  possible  alloca¬ 
tions  in  order  to  find  the  exact  optimal  solution.  This  scenario  was  considered  previously  for 
the  deterministic  case  in  the  limit  of  large  population  sizes  ( N 1  =  100,000  and  N2  =  200,000, 
using  our  notation)  by  Keeling  and  Shattock  [2],  who  found  that  for  a  wide  range  of  values  of 
the  reproductive  number  r0,  the  optimal  solution  as  a  function  of  increasing  V  was  governed  by 
the  ability  to  induce  herd  immunity  in  the  smaller  population  (small  V),  then  in  the  larger  pop¬ 
ulation  (intermediate  V),  and  finally  in  both  (large  V).  The  colormap  of  Figs.  4A  and  4B  illus¬ 
trates  the  average  epidemic  size  (E)  corresponding  to  a  particular  allocation  of  vaccine, 
quantified  in  the  figure  by  the  amount  of  vaccine  in  the  smaller  population  V1.  The  optimal 
vaccine  allocation  is  illustrated  by  the  black  line,  along  which  (E)  is  minimized.  The  range  in 
y1  is  limited  by  the  constraints  V1  <  V  —  /p  =  V  —  1,  and  Vt  >  V  -  999,  i.e.  neither  popula¬ 
tion  receives  more  vaccine  than  the  number  of  initial  susceptible  individuals  in  that  population. 
This  results  in  the  limiting  diagonals  in  the  colormap  of  Figs.  4A  and  4B. 

Switching 

In  the  deterministic  model,  Keeling  and  Shattock  [2]  found  that  the  optimal  solution  exhibited 
“switching”  behavior,  in  which  the  optimal  vaccine  allocation  makes  a  significant,  discontinu¬ 
ous,  change  when  the  total  amount  of  vaccine  V  exceeds  a  threshold  size.  When  the  amount  of 
vaccine  V  is  below  this  threshold  size,  the  majority  of  the  vaccine  is  optimally  allocated  to  the 
smaller  population.  When  the  amount  of  vaccine  V  is  above  this  size,  all  of  it  is  optimally  allo¬ 
cated  to  the  larger  population.  This  behavior  persists  for  a  wide  range  of  reproductive  numbers 
r0  >  1  in  the  deterministic  SIR  model. 

The  stochastic  optimal  solution  exhibits  switching  behavior  only  for  larger  values 
of  r0  (r0  a  2.9,  for  the  other,  fixed  parameters  considered  here).  This  is  demonstrated  for  r0  =  5 
in  Fig.  4A.  Switching  occurs  first  at  V  =  474  and  then  again  at  V  =  780.  The  first  switching 
point,  above  which  all  vaccine  is  optimally  allocated  to  the  larger  population,  is  present  in  both 
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the  stochastic  and  deterministic  models,  although  the  switching  point  of  the  stochastic  model 
occurs  at  a  smaller  amount  of  vaccine  V.  The  second  switching  point  is  absent  in  the  determin¬ 
istic  model.  The  presence  of  the  second  switch  in  the  stochastic  model  is  explained  in  terms  of 
the  relative  heights  of  the  peaks  of  the  gain  curves  in  the  next  subsection. 

For  intermediate  values  of  r0  (1  <  r0  Si  2.9)  in  the  stochastic  model,  there  is  no  switching  be¬ 
havior,  which  is  in  contrast  to  the  results  of  the  deterministic  model.  It  is  instead  optimal  to  dis¬ 
tribute  any  given  total  amount  of  vaccine  V  approximately  in  proportion  to  the  sizes  of  the 
populations  themselves.  This  is  shown  for  r0  =  2  in  Fig.  4B.  The  continuous  transition  between 
large  r0  values  where  switching  does  take  place  and  small  r0  values  where  it  does  not  is  illustrat¬ 
ed  in  Fig.  4C.  As  r0  is  decreased,  the  region  between  the  two  instances  of  switching  behavior, 
where  all  vaccine  is  taken  out  of  the  smaller  population  and  V1  =  0,  becomes  narrower  and  dis¬ 
appears  completely  between  r0  =  2  and  r0  =  3  (at  r0  «  2.9). 

The  conclusion  for  the  stochastic  optimal  solution  is  that  for  intermediate  values 
of  r0  ( 1  <  r0  <  2.9),  the  optimal  solution  is  to  approximately  distribute  vaccine  in  proportion 
to  population  size.  For  large  r0  (r0  a  2.9),  two  switches  take  place.  In  the  deterministic  optimal 
solution,  a  single  switch  takes  place  for  all  values  of  r0  >  1.  Note  that  the  stochastic  gain  curves 
exhibit  peaks  only  when  r0  S  2.5,  while  the  deterministic  gain  curves  always  exhibit  peaks.  The 
relationship  between  switching  and  the  presence  of  peaks  in  the  gain  curve  will  be  discussed  in 
the  following  subsection. 

Understanding  the  Optimal  Stochastic  Solution 

Next  we  examine  the  optimal  stochastic  solutions  of  Fig.  4A  and  4B  in  closer  detail.  Figs.  5A 
and  5B  show  the  optimal  stochastic  solution  as  a  solid  line  and  the  optimal  deterministic  solu¬ 
tion  as  a  dashed  line,  with  both  solutions  represented  by  the  fraction  of  the  amount  of  vaccine 
V1/ V  given  to  the  smaller  population.  We  seek  to  characterize  the  different  strategies  employed 
by  the  optimal  solution  in  the  different  resource  regimes,  and  also  to  quantify  why  the  alloca¬ 
tion  transitions  abruptly  from  one  population  to  another.  More  broadly,  we  explain  why  the 
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Fig  5.  Optimal  Protocols  Represented  as  a  Fraction  of  Total  Vaccine.  The  solid  blue  line  shows  the  optimal  fraction  of  the  available  vaccine  allocated  to 
the  smaller  population  t/Vl/in  order  to  minimize  the  average  stochastic  epidemic  size  (E).  The  dashed  blue  line  shows  a  different  optimal  allocation  that 
minimizes  the  deterministic  epidemic  size  E.  The  right  vertical  scale  applies  to  these  two  measurements.  The  difference  in  the  resultant  average  epidemic 
size  {£}  between  the  two  protocols  when  applied  to  the  stochastic  model  is  plotted  in  gold  (measured  by  the  left  vertical  scale).  The  results  are  obtained  for  r0 
=  5  in  Figure  A,  and  r0  =  2  in  Figure  B,  and  populations  of  sizes  N 1  =  500  individuals  and  N 2  =  1 000  individuals,  both  initialized  with  a  single  infected  individual 
I'  =  I2  =  i 
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optimal  stochastic  solution,  which  minimizes  the  average  epidemic  size  (E) ,  differs  from  the 
optimal  deterministic  solution,  which  minimizes  the  characteristic  epidemic  size  E. 

Much  of  our  insight  comes  from  analysis  of  the  gain  curves,  which  are  shown  in  Fig.  3A  and 
3B  for  r0  =  5  and  r0  =  2.  By  Equation  6,  the  area  under  each  gain  curve  G(V)  up  to  a  particular 
value  of  V  is  the  decrease  in  the  epidemic  size  due  to  that  amount  of  vaccine  V.We  begin  with 
the  case  when  r0  =  5,  which  is  representative  of  the  large  r0  regime  (r0  a  2.9)  where  switching 
does  take  place. 

Fig.  5  A  shows  that  with  small  amounts  of  vaccine  V,  all  of  the  vaccine  is  optimally  allocated 
to  the  smaller  population.  In  the  deterministic  case,  this  strategy  persists  for  larger  amounts  of 
vaccine  V,  unless  there  is  enough  vaccine  for  the  small  population  to  achieve  herd  immunity,  at 
V  =  400.  For  a  range  of  V  greater  than  V  =  400,  herd  immunity  is  preserved  in  the  small  popu¬ 
lation,  and  the  remaining  vaccine  is  optimally  allocated  to  the  large  population.  For  the  sto¬ 
chastic  case,  vaccine  allocation  to  the  larger  population  begins  for  a  smaller  V,V=  324, 
above  which  the  optimal  solution  is  to  maintain  V]  =  324  and  devote  the  remainder  of  the  vac¬ 
cine  to  V2. 

This  difference  in  strategy  can  be  attributed  to  the  fact  that  in  the  stochastic  model,  the  gain 
G  begins  to  decline  well  before  the  onset  of  effective  herd  immunity.  This  is  evident  in  both  the 
N 1  and  N2  solid  curves  of  Fig.  3  A.  In  contrast,  the  gain  in  the  deterministic  model  peaks  very 
close  to  herd  immunity  at  V  =  400,  as  the  dashed  curves  of  Fig.  3  A  show.  Compared  to  the  de¬ 
terministic  model,  one  can  attribute  this  earlier  decline  in  the  average  epidemic  size  (E)  as 
being  due  to  the  probability  of  a  terminal  infection,  which  significantly  lowers  the  average  (E) . 
More  quantitatively,  in  Fig.  3A,  V  =  324  is  the  point  at  which  the  stochastic  curve  for  N1  crosses 
the  initial  value  of  curve  N2. 

As  V  increases  further  in  Fig.  5A,  there  is  a  sharp  transition,  indicating  that  if  more  vaccine 
exists  than  V  =  474  in  the  stochastic  model  or  V  =  657  in  the  deterministic  model,  all  vaccine 
should  optimally  be  allocated  to  the  large  population.  This  is  the  first  switch  noted  earlier  that 
occurs  in  both  models.  As  with  the  earlier  transition,  the  first  switch  takes  place  at  a  smaller 
amount  of  vaccine  V  in  the  stochastic  model  than  in  the  deterministic  model.  This  sacrifices 
herd  immunity  that  could  have  been  achieved  in  the  small  population,  in  favor  of  relatively 
larger  gains  in  protection  that  can  be  achieved  with  this  level  of  vaccine  in  the  large  population. 
Quantitatively  it  is  clear  from  Fig.  3A  that  beyond  a  certain  amount  of  vaccine,  the  stochastic 
gain  curve  G{V)  begins  to  decline  for  N2  while  the  stochastic  curve  for  N 1  is  still  relatively  large 
and  constant.  Thus  around  this  level  of  vaccine,  all  the  available  vaccine  should  optimally  be 
switched  into  the  larger  population.  This  same  behavior  is  observed  for  the  deterministic  gain 
curves  G(  V)  for  correspondingly  larger  values  of  V. 

Complete  resource  allocation  to  the  large  population  continues  until  the  large  population 
achieves  herd  immunity,  at  which  point  a  fraction  of  the  vaccine  is  allocated  to  the  smaller  pop¬ 
ulation.  For  the  deterministic  case,  the  optimal  solution  retains  herd  immunity  for  the  large 
population,  and  increasingly  allocates  resources  to  the  small  population,  until  both  populations 
achieve  herd  immunity.  After  that  point,  the  optimal  solution  plateaus.  For  the  deterministic 
model,  the  epidemic  never  progresses  (1(f)  <  I0).  Because  there  is  no  cost  for  vaccination,  re¬ 
maining  resources  are  allocated  based  solely  on  the  relative  population  sizes  (i.e.  1/3  for  the 
small  population  and  2/3  for  the  large  population).  For  the  deterministic  model,  this  corre¬ 
sponds  to  a  situation  with  excess  vaccine,  since  both  populations  are  fully  protected  once  each 
has  sufficient  resources  to  insure  herd  immunity. 

For  the  stochastic  model,  once  there  is  sufficient  vaccine  to  induce  effective  herd  immunity 
in  the  large  population,  at  around  V  =  660,  vaccine  is  once  again  allocated  to  the  small  popula¬ 
tion.  However,  unlike  the  deterministic  case,  for  the  stochastic  model,  there  is  a  second  abrupt 
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shift  in  resources  around  V  =  780,  resulting  in  a  cusp  in  the  optimal  V1/V,  with  the  optimal  so¬ 
lution  approaching  the  final  population  based  plateau  value  V1/  V  =1/3  from  above. 

This  is  due  to  the  fact  that  in  the  stochastic  model,  for  large  r0,  the  smaller  population  N1 
has  a  greater  peak  in  gain.  Thus  if  there  is  enough  vaccine  available,  there  is  a  benefit  to  remov¬ 
ing  some  vaccine  from  the  large  population  in  order  to  take  advantage  of  the  higher  gain  in  the 
smaller  population.  This  second  switch  does  not  occur  in  the  deterministic  model  because  the 
opposite  is  true,  the  peak  of  deterministic  curve  N2  for  the  larger  population  is  always  higher 
than  the  peak  of  deterministic  curve  N1  for  the  smaller  population  in  Fig.  3A. 

For  the  stochastic  model,  the  gain  curves  of  Fig.  3B  can  also  be  used  to  explain  the  absence 
of  switching  behavior  for  r0  =  2,  which  is  generally  observed  for  lower  values  of  r0  (1  <  r0  Si 
2.9).  A  significant  difference  in  this  case  is  that  the  gain  decreases  continuously.  Due  to  the  ab¬ 
sence  of  peaks  in  the  gain  curve,  the  second  switch  observed  for  the  large  r0  stochastic  model, 
does  not  occur  for  small  values  of  r0.  The  absence  of  the  first  switch  is  more  subtle  and  depends 
on  more  than  just  the  presence  of  a  peak  which  exists  when  r0  S  2.5  as  discussed  previously. 
For  the  first  switch  to  occur,  the  peak  of  the  curve  G(V)  must  be  large  enough  to  offset  the  de¬ 
clines  in  the  gain  that  first  population  N 1  exhibits.  Hence  the  first  switch  takes  place  for  a  more 
restrictive  set  of  r0,  and  only  when  the  peak  in  the  gain  is  sufficiently  large  (r0  &  2.9). 

In  summary,  the  switching  behavior  of  the  optimal  vaccination  allocation  are  due,  firstly,  to 
the  presence  of  peaks  in  the  gain  curves,  and  secondly,  due  to  the  relative  heights  of  these 
peaks.  This  explains  why  in  the  stochastic  model,  switching  occurs  only  for  large  values  of  r0, 
while  in  the  deterministic  model,  it  occurs  for  all  values  of  r0  >  1 .  Fundamentally,  this  differ¬ 
ence  arises  from  the  bimodal  nature  of  the  epidemic  size  distribution  P{E). 


Conditioning  on  the  Large-Scale  Epidemic 

In  the  previous  sections  we  identified  and  analyzed  the  optimal  vaccine  allocations,  obtained 
by  minimizing  the  average  epidemic  size  (E),  when  both  populations  are  initialized  with  a  sin¬ 
gle  infected  individual  I*  =  I2  =  1.  For  such  small  initial  infection  exposures,  the  average  epi¬ 
demic  size  (E)  has  significant  contributions  from  the  terminal  infection  and  the  large-scale 
epidemic  of  both  populations. 

Alternatively,  a  policymaker  may  be  interested  in  a  more  cautious  approach  which  ignores 
contributions  from  terminal  infections  in  the  optimization,  focusing  instead  on  optimal  alloca¬ 
tion  assuming  large-scale  epidemics  are  likely  to  develop.  This  scenario  is  achieved  in  the  sto¬ 
chastic  model  by  initializing  the  system  with  a  relatively  large  number  of  initial  infected 
individuals  7*  and  I2.  Results  are  shown  below  for  two  populations  with  N 1  =  500  and  N2  = 

1000  individuals  and  7<5  =  =  70  initially  infected  individuals. 

Both  Figs.  6A  and  6B  illustrate  that  as  the  initial  number  of  infected  individuals  70  increases, 
the  stochastic  optimal  distribution  curves  (solid  lines)  are  increasingly  similar  to  the  determin¬ 
istic  optimal  protocol  (dotted  line),  both  qualitatively  and  quantitatively.  In  the  r0  =  5  case,  the 
second  vaccine  switch  at  V  =  780,  characteristic  of  the  stochastic  optimum,  disappears  between 
70  =  1  and  70  =  2.  Similarly  for  r0  =  2,  the  optimal  distribution  changes  incrementally  from  an 
approximately  proportional  distribution  when  70  =  1  (blue  line),  to  a  distribution  with  a  single 
switch  when  70  =  5  (green  line)  that  has  features  similar  to  those  of  the  optimal  distribution 
curve  for  the  deterministic  model. 

For  the  case  of  two  populations,  Nl  and  N2,  we  define  the  overall  gain  G  to  be  the  magnitude 
of  the  incremental  decrease  in  the  average  combined  epidemic  size  (E)  per  incremental  increase 
in  the  total  amount  of  vaccine  V.  Figs.  6C  and  6D  illustrate  the  overall  gain  curves  associated 
with  each  of  the  optimal  solutions  illustrated  in  Figs.  6A  and  6B,  respectively.  Note  that  the  over¬ 
all  gain  for  the  two  population  optimization  problem  generalizes  the  notion  of  gain  for  an 
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Fig  6.  Variation  of  the  Optimal  Allocation  with  Increasing  Likelihood  of  a  Large-Scale  Epidemic.  The  solid  lines  in  Figures  A  and  B  illustrate  the  optimal 
fraction  of  the  available  vaccine  allocated  to  the  smaller  population  V'lV  in  order  to  minimize  the  average  stochastic  epidemic  size  <£)  for  varying  values  of 
l0  =  l'0  =  Pa.  The  dashed  black  line  shows  the  optimal  allocation  that  minimizes  the  deterministic  epidemic  size  E.  Figures  C  and  D  illustrate  the  corresponding 
overall  gain  G  as  a  function  of  the  total  vaccine  V  for  each  of  the  curves  in  Figures  A  and  B.  The  results  are  obtained  for  r0  =  5  in  Figures  A  and  C,  and/'0  =  2  in 
Figures  B  and  D,  and  populations  of  sizes  A/1  =  500  individuals  and  N2  =  1 000  individuals.  As  the  initial  seed  population  of  infected  individuals  is  increased, 
from  l'0  =  P0  =  1  (as  in  Fig.  4),  to  l'0  =  t20  =  5,  both  populations  are  increasingly  likely  to  experience  a  large-scale  epidemic.  Specifically,  in  the  stochastic 
model,  the  probability  of  having  large-scale  epidemics  in  both  populations  for  /0  =  1 , 2,  5  is,  respectively,  0.6392  (blue),  0.921 0  (red),  and  0.9993  (green)  for 
r0  =  5  and  0.2468,  0.5559,  and  0.9334  for  r0  =  2.  While  both  the  optimal  allocation  and  gain  in  the  stochastic  model  becomes  more  similar  to  the  deterministic 
limit  with  increasing  /0,  it  does  not  converge  to  the  deterministic  limit  for  finite  r0  due  to  the  non-negligible  width  of  the  distribution  of  large-scale  epidemics  that 
is  observed  in  the  stochastic  case. 
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individual  population  that  was  introduced  in  Equation  12  and  Fig.  3.  For  a  single  population, 
gain  varies  smoothly  with  increasing  V.  In  contrast,  for  two  populations,  the  overall  gain  curve 
incorporates  portions  of  the  gain  curves  for  the  individual  populations,  and  may  exhibit  sharp 
kinks  and  discontinuous  jumps,  reflecting  changes,  such  as  switching,  which  occur  at  transition 
points  where  the  allocation  changes  discontinuously  from  one  population  to  the  other. 

Comparing  the  overall  gain  curves  in  Figs.  6C  an  6D  with  the  corresponding  optimal  proto¬ 
cols  in  Figs.  6A  and  6B,  we  see  that  switching  (i.e.  discontinuous  jumps  in  the  optimal  protocol) 
coincides  with  a  discontinuous  increase  in  gain.  Kinks  (slope  discontinuities)  in  the  optimal 
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solution  (i.e.  points  where  the  protocol  shifts  from  complete  allocation  to  one  population,  grad¬ 
ually  increasing  the  allocation  of  the  other  population)  coincide  with  kinks  at  local  minima  of 
the  overall  gain.  Peaks  in  the  overall  gain  occur  at  intermediate  points,  rather  than  turning 
points,  in  the  protocol.  As  the  initial  number  of  infected  individuals  I0  increases  in  each  of  the 
populations,  the  amplitude  of  the  peaks  in  the  overall  gain  curves  become  increasingly  sharp, 
becoming  more  similar,  yet  not  identical  to,  the  corresponding  overall  gain  of  the  deterministic 
optimal  solution. 

For  finite  values  of  r0  such  as  those  considered  here,  the  optimal  distribution  curves  also  ap¬ 
proach  but  never  fully  converge  to  the  deterministic  optimal  solution.  Comparing  the  I0  =  5 
(green  lines)  and  deterministic  curves  (dotted  lines),  appreciable  differences  remain  apparent, 
even  though  the  probability  of  a  large-scale  epidemic  exceeds  99.9%  for  r0  =  5,  and  93.3%  for  r0 
=  2.  Even  after  conditioning  on  large-scale  epidemics,  the  non-zero  width  of  the  probability 
distribution  for  large-scale  epidemics  that  is  present  in  the  stochastic  models,  and  absent  in  the 
deterministic  case,  factors  nontrivially  into  the  optimization  problem.  This  is  the  primary  fac¬ 
tor  contributing  to  the  difference  between  the  stochastic  and  deterministic  optimal  solutions. 

We  tested  even  larger  values  of  70  (e.g.  I0  =  10,  not  shown)  to  verify  the  non-convergence  be¬ 
tween  the  stochastic  and  deterministic  protocols  in  the  large  I0  limit.  However,  the  system 
states  generated  by  increasing  I0,  leaving  S0  +  I0  =  N  fixed,  become  unrealistic  if  I0  is  taken  to  be 
too  large.  In  fact,  the  optimal  allocation  curves  for  I0  =  10  lie  farther  from  the  deterministic 
curves  (dotted  lines)  than  the  I0  =  5  curves  (green  lines).  This  is  due  to  the  fact  that  there  is 
only  a  very  small  probability  of  a  system  initialized  with  I0  =  1  and  S0  =  999  individuals,  as  in 
the  deterministic  model,  to  later  transition  into  a  state  with  I0  =  10  and  So  =  990  individuals 
without  the  recovery  of  any  infected  individuals  at  intermediate  times.  Thus,  while  increasing 
the  number  of  initially  infected  individuals  I0  brings  the  system  closer  to  the  deterministic  limit 
at  first,  initializing  the  system  with  too  many  initial  infected  I0  creates  an  unrealistic  scenario. 

The  Range  of  Outcomes 

The  differences  in  the  optimal  vaccination  protocols  between  the  stochastic  and  deterministic 
models  can  lead  to  substantial  differences  in  the  observed  outcomes.  The  optimal  protocols  for 
the  stochastic  and  deterministic  models  coincide  for  small  quantities  of  vaccine,  where  in  both 
cases  it  is  optimal  to  allocate  all  vaccines  to  the  smaller  population.  The  optimal  solutions  also 
coincide  in  the  limit  of  large  quantities  of  vaccine,  where  it  is  optimal  to  allocate  vaccines  in 
proportion  to  the  population  size.  The  range  of  possible  outcomes  for  different  vaccination 
strategies  indicates  that  optimization  is  most  important  when  intermediate  amounts  of  vaccine 
are  available.  One  way  of  understanding  the  potential  impact  associated  with  the  optimal  sto¬ 
chastic  and  deterministic  protocols  is  by  comparing  their  projected  outcomes  when  applied  to 
the  presumably  more  realistic  stochastic  SIR  model.  In  this  scenario,  there  is  considerable  dif¬ 
ference  between  the  best  and  worst  possible  outcomes  and  a  significant  but  smaller  difference 
between  the  stochastic  and  deterministic  optimal  solutions. 

The  dashed  gold  lines  of  Figs.  5  A  and  5B  illustrate  the  difference  in  the  resultant  average  ep¬ 
idemic  size  (E)  between  the  stochastic  and  deterministic  optimal  protocols.  Both  protocols  are 
the  same  in  resource  rich  and  resource  poor  regimes,  and  hence  yield  identical  results.  Figs.  7  A 
and  7B  illustrate  (E)  for  both  the  stochastic  and  deterministic  optimal  protocols  as  well  as  the 
worst  case  allocation. 

We  define  the  “worst  case”  protocol  as  that  which  maximizes  (E)  within  the  range  of  al¬ 
lowed  allocations  illustrated  in  Fig.  4A  and  4B.  Together  the  stochastic  optimal  solution  and 
worst  case  allocation  define  the  possible  range  of  (E)  at  a  given  value  of  V.  Fig.  7  shows  that  the 
difference  between  the  outcome  of  the  worst  case  protocol,  and  either  the  optimal  stochastic 
and  deterministic  protocols,  is  substantially  larger  than  the  difference  between  the  stochastic 
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Fig  7.  Comparison  of  Average  Epidemic  Size  for  Optimal,  Deterministic  Optimal,  and  Worst  Case  solutions.  These  figures  illustrate  the  average 
epidemic  size  (E)  forthree  different  protocols:  the  stochastic  optimum  (minimizes  (£)),  the  deterministic  optimum,  and  the  worst  case  scenario  (maximizes 
(£)).  The  three  different  colors  correspond  to  different  number  of  initial  infected  individuals  /0.  Results  are  obtained  for  two  populations  with  size  A/1  =  500 
individuals  and  N2  =  1000  individuals.  Results  are  shown  forr0  =  5  in  Figure  A,  and/-0  =  2  in  Figure  B.  Other  than  the  value  of  r0,  the  color  and  line  style 
legends  in  Figure  B  apply  to  both  graphs. 
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and  deterministic  cases.  This  is  particularly  pronounced  for  smaller  values  of  r0,  i.e.  r0  =  2.  The 
worst  case  protocol  would  involve  continuing  to  place  vaccine  in  a  population  even  after  it  is 
near  or  has  reached  herd  immunity.  This  is  represented  by  the  plateaus  where  the  average  epi¬ 
demic  size  (E)  is  not  significantly  lowered  by  further  vaccinating  members  of  the  population. 
Deterministically,  this  is  evident  from  the  fact  that  d  l/dt  <  0  as  soon  as  the  herd  immunity 
threshold  has  been  reached.  The  deterministic  herd  immunity  threshold,  serves  as  an  approxi¬ 
mate  guide  for  when  to  stop  vaccinating  even  in  the  stochastic  case. 

The  differences  between  the  stochastic  and  deterministic  protocols  have  a  complex  r0  and  I0 
dependence.  The  effect  on  the  difference  in  (E)  between  the  stochastic  and  deterministic  mod¬ 
els  that  is  caused  by  increasing  I0  is  different  for  small  compared  to  large  reproductive  numbers 
r0.  With  a  small  reproductive  number,  e.g.,  r0  =  2,  the  difference  in  the  average  epidemic  size 
between  the  stochastic  and  deterministic  optimal  protocols  is  largest  at  an  intermediate  value 
of  I0,  Ig  =  1^  =  2  for  the  case  illustrated  in  Fig.  7B.  In  contrast  for  large  reproductive  number, 
e.g.  r0  =  5,  Fig.  7  A  illustrates  that  the  difference  in  the  average  epidemic  size  between  the  sto¬ 
chastic  and  deterministic  optimal  protocols  is  maximized  for  Ig  =  Ig  =  1  and  decreases  steadi¬ 
ly  as  I0  is  increased. 

Comparison  with  Proportional  Distribution 

In  the  previous  section  we  compared  average  epidemic  sizes  (E)  obtained  in  the  stochastic 
model  when  the  stochastic,  deterministic,  and  worst  case  outcome  protocols  are  applied.  An¬ 
other  relevant  protocol  for  comparison  is  one  where  vaccine  is  always  distributed  in  a  manner 
that  is  proportional  to  the  population  sizes,  i.e.  V1  =  {N^Ni ;  V  and  V2  =  N^N,  V.  Politically,  a 
proportional  distribution  of  vaccine  is  expected  to  be  much  easier  to  implement  with  the  pub¬ 
lic,  compared  to  a  protocol  that  allocates  vaccine  to  one  population  at  the  expense  of  another, 
even  if  the  expected  epidemic  size  is  reduced  in  the  skewed  distribution. 

Fig.  8  illustrates  the  increase  in  the  average  epidemic  size  (E)  for  a  proportional  distribution 
of  vaccine  in  the  stochastic  model,  compared  to  the  optimal  stochastic  solution  for  different 
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Fig  8.  Comparison  of  Proportional  and  Stochastic  Optimal  Protocols.  The  difference  in  average  epidemic  size  (E)  obtained  in  the  stochastic  model,  for  a 
distribution  that  is  directly  proportional  to  population  size  (A/1  =  500  and  N2=  1000),  compared  to  the  optimal  stochastic  protocol.  Results  are  shown  forr0  =  5 
(Figure  A),  and  r0  =  2  (Figure  B),  and  different  values  of  of  l0  (l0  =  1 , 2,  and  5).  The  differences  are  most  significant  for  large  values  of  r0  and  /0,  where  the 
behavior  of  the  stochastic  model  is  most  similar  to  the  deterministic  case.  At  smaller  values  of  r0  and  /0  the  distinction  is  much  less  pronounced.  The 
difference  approaches  zero  for  large  amounts  of  vaccine,  after  effective  herd  immunity  is  reached  in  both  populations,  and  the  optimal  protocol  approaches 
proportional  distribution.  Note  the  difference  in  scale  on  the  vertical  axes  for  Fig.  8A  and  8B. 
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values  of  r0  and  I0.  Alternatively,  the  curves  can  be  interpreted  as  the  average  reduction  in  the 
epidemic  size  when  the  distribution  changes  from  a  proportional  to  the  stochastic  optimal  dis¬ 
tribution.  The  most  significant  differences  occur  for  large  values  of  r0  and  70  (i.e.  r0  =  5),  which 
is  the  regime  in  which  the  stochastic  model  is  most  similar  to  the  deterministic  model.  In  con¬ 
trast,  when  r0  and  I0  are  relatively  small  (i.e.  r0  =  2),  the  difference  between  the  stochastic  and 
proportional  distributions  is  much  less  significant.  Additionally,  for  all  values  of  the  parame¬ 
ters,  when  there  is  sufficient  vaccine  available  to  ensure  effective  herd  immunity  in  both  popu¬ 
lations,  the  optimal  protocol  approaches  the  proportional  distribution,  and  the  difference  in 
the  average  epidemic  size  for  the  two  cases  approaches  zero. 

Population  Size  Variations 

Thus  far  we  have  considered  the  case  of  two  noninteracting  populations,  where  one  population 
is  twice  as  large  as  the  other,  i.e.,  Nl  =  500  and  N2  =  1000.  In  this  section,  we  show  that  the  qual¬ 
itative  characteristics  of  the  optimal  stochastic  and  deterministic  protocols  are  robust  to  varia¬ 
tions  in  the  ratio  of  the  population  sizes  A/V/V2.  For  this  investigation,  the  value  of  N2  =  1000  is 
held  constant,  and  we  compute  optimal  stochastic  protocols  by  minimizing  (E)  for  N1  =  250, 
500,  750,  and  1000. 

Fig.  9  illustrates  the  optimal  stochastic  (Figs.  9A  and  9B)  and  deterministic  (9C  and  9D) 
protocols  for  r0  =  5  (9 A  and  9C)  and  r0  =  2  (9B  and  9D)  for  varying  values  of  N1.  The  same 
qualitative  behavior  observed  in  Figs.  4  and  5  for  N 1  =  500  is  illustrated  here  for  different  ratios 
of  the  population  sizes.  For  r0  =  5  the  stochastic  and  deterministic  models  exhibit  switching  in 
a  manner  that  is  qualitatively  similar  to  our  previous  results  (two  switches  in  the  stochastic 
case,  and  one  in  the  deterministic  case).  And  as  before,  when  r0  =  2,  switching  is  absent  in  the 
stochastic  model,  but  preserved,  with  one  switch,  in  the  deterministic  case. 

The  only  exception  arises  when  N 1  =  1000,  i.e.  N]  /N2  =  1.  In  this  case,  there  is  no  broken 
symmetry  in  the  population  sizes.  In  the  deterministic  model,  for  both  r0  =  5  and  r0  =  2,  the  op¬ 
timal  protocol  allocates  vaccine  to  one  population  (either  population  can  be  chosen),  until  herd 
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Fig  9.  Impact  of  Variations  of  the  Population  Size  Ratio  on  the  Optimal  Protocol.  Figures  A  and  B  illustrate  results  for  the  stochastic  model  with  r0  =  5 
and  r0  =  2,  respectively.  Figures  C  and  D  illustrate  corresponding  results  for  the  deterministic  model.  Each  curve  shows  the  optimal  fraction  of  the  available 
vaccine  allocated  to  the  smaller  population  I/1/!/  in  order  to  minimize  the  average  epidemic  size  (E)  of  the  population  as  a  whole.  In  all  cases,  the  number  of 
individuals  in  the  second  population  N2  =  1 000  is  kept  constant.  Different  curves  illustrate  results  as  the  number  of  individuals  in  the  first  population  is  varied: 
A/1  =  250,  500,  750,  and  1 000.  The  qualitative  results  obtained  previously  for  the  stochastic  and  deterministic  protocols  in  Figs.  4  and  5  are  preserved  for 
changes  in  the  ratio  of  population  sizes,  except  in  the  case  of  equal  population  sizes  N1  IN2  =  1 ,  in  which  case  there  is  no  broken  symmetry  in  the  population 
sizes.  All  populations  are  initialized  with  a  single  infected  individual  I]  =  /§  =  1 . 
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immunity  is  reached.  At  higher  levels  of  vaccine,  the  deterministic  protocol  allocates  the  re¬ 
maining  vaccine  to  the  other  population  until  it  also  achieves  herd  immunity.  Subsequently, 
vaccine  is  divided  proportionally.  For  the  stochastic  model,  when  r0  =  5  and  N1  =  1000 
(Fig.  9A),  a  protocol  that  is  similar  to  the  deterministic  case  is  optimal,  but  involves  effective 
herd  immunity,  and  transitions  earlier  and  more  sharply  to  a  proportional  distribution.  How¬ 
ever,  when  r0  =  2  in  the  stochastic  model  (Fig.  9B),  the  optimal  solution  transitions  into  exactly 
the  proportional  distribution  beginning  at  a  very  small  value  of  V. 

Alternate  Cost  Functions 

So  far,  we  have  defined  the  optimal  allocation  as  that  which  minimizes  the  average  epidemic 
size  (E),  a  quantity  that  contains  contributions  from  both  terminal  infections  and  large-scale 
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epidemics,  but  is  not  necessarily  representative  of  any  specific  epidemic  size  that  is  likely  to  be 
observed  because  of  the  gap  in  the  size  distribution  P{E)  (Fig.  1A).  Choosing  to  minimize  the 
deterministic  result,  which  is  the  same  as  the  average  large-scale  epidemic  size,  might  potential¬ 
ly  be  viewed  as  a  conservative  approach  that  safeguards  against  the  case  in  which  both  popula¬ 
tions  will  experience  large-scale  epidemics. 

Other  criteria  for  optimization  may  be  considered  within  this  framework.  For  example,  one 
possibility  would  be  to  maximize  the  likelihood  of  having  no  (or  very  few)  infections.  We  com¬ 
puted  the  distribution  of  vaccine  that  would  maximize  the  probability  of  having  no  further  in¬ 
fection  beyond  the  initial  infectious  seed.  In  this  case,  the  optimal  protocol  allocates  all  vaccine 
to  the  smaller  population  until  every  individual  is  vaccinated,  only  allocating  vaccine  to  the 
larger  population  when  V  exceeds  N1. 

A  cost  function  that  may  be  of  particular  interest  to  policymakers  is  one  that  illustrates  how, 
for  a  given  quantity  of  vaccine  V,  to  allocate  vaccine  so  as  to  minimize  the  probability  of  having 
an  epidemic  greater  than  a  particular  size.  To  address  this,  in  the  same  scenario  of  two  non-in¬ 
teracting  populations  (e.g.  two  well-separated  cities)  with  N 1  =  500  individuals  and  N2  =  1000 
individuals  and  both  populations  initialized  with  a  single  initial  infected  7*  =  I2  =  1,  here  we 
alternatively  consider  the  probability  that  the  epidemic  is  below  some  particular  threshold  tol¬ 
erance  size  Emax. 

A  policymaker  may  be  interested  in  how  much  vaccine  V  would  be  necessary  and  how  it 
must  be  allocated  between  two  populations  in  order  to  keep  the  total  epidemic  below  some  size 
£max.  We  compute  the  best  achievable  probability  of  having  an  epidemic  below  a  given  size 
£max  given  a  total  amount  of  vaccine  V.  The  results  are  shown  in  Fig.  10  A. 

The  sharp  color  contrast  of  the  diagonal  bands  in  Fig.  10A  are  associated  with  step-like 
changes  in  probability,  arising  from  the  bimodal  nature  of  the  epidemic  size  distributions  P(E). 
Because  there  is  very  little  probability  for  an  event  in  the  size  range  between  the  large-scale  epi¬ 
demic  and  the  terminal  infection  peaks,  when  the  threshold  Emax  passes  through  the  large-scale 
epidemic  size  (which  depends  on  the  vaccine  allocation)  in  the  small  population,  the  large  pop¬ 
ulation,  or  the  sum  of  the  two,  nearly  discrete  steps  in  probability  are  observed. 

The  allocation  that  maximizes  this  probability  is  shown  in  Fig.  10B,  and  is  a  function  of 
both  the  amount  of  vaccine  V  and  also  £max.  Unlike  our  previous  optimization  based  on  ex¬ 
pected  size  (where  the  corresponding  plot  depends  only  on  V),  here  the  solution  is  extremely 
complex,  switching  discontinuously  and  frequently  depending  on  both  V  and  £max,  as  indicat¬ 
ed  by  sharp  grey  scale  contrasts  reflecting  boundaries  between  high  and  low  allocations  to  the 
small  population.  In  the  resource  poor  regime  (small  V,  corresponding  to  the  lower  horizontal 
boundary  of  the  color  plot)  the  solution  switches  from  full  allocation  to  the  small  population, 
to  full  allocation  to  the  large  population,  back  to  full  allocation  to  the  small  population.  The 
lower  left  white  triangle  in  Fig.  10B  corresponds  to  the  situation  with  few  resources,  and  mini¬ 
mal  tolerance  for  the  epidemic  size.  As  in  the  previous  stochastic  and  deterministic  solutions 
aimed  at  minimizing  the  average  epidemic  size,  here  the  optimal  solution  allocates  all  resources 
to  the  smaller  population.  In  the  £max  dependent  resource  rich  regime,  corresponding  to  points 
above  the  highest  diagonal,  the  maximum  achievable  probability  in  Fig.  10A  is  near  unity,  and 
the  optimal  allocation  simplifies  to  depend  only  on  V  (corresponding  to  horizontal  bands  in 
Fig.  10B).  However,  in  intermediate  cases,  where  tradeoffs  are  most  critical,  the  structure  of  the 
resulting  solution  is  much  too  subtle  to  be  realistically  implemented  for  real  populations  given 
a  limited  amount  of  vaccine  V. 

For  comparison,  we  evaluate  the  corresponding  probabilities  based  on  our  previous  stochas¬ 
tic  and  deterministic  optimal  protocols.  While  both  solutions  are  suboptimal  for  this  alterna¬ 
tive  criterion,  the  stochastic  solution  comes  close  to  the  optimal  case.  Fig.  10C  shows  this  result 
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Fig  1 0.  Optimizing  the  probability  of  having  an  epidemic  less  than  a  given  size.  Figure  A  shows  in  color,  the  optimal  (largest)  probability  of  having  an 
epidemic  less  than  some  given  size  (x-axis),  given  some  amount  of  vaccine  (y-axis).  Figure  B  shows  the  fraction  of  vaccine  I/1  IV  in  the  smaller  population  that 
corresponds  to  optimal  probability  shown  in  Figure  A.  Figures  C  and  D  show  the  same  probabilities  for,  respectively,  the  deterministic  and  the  stochastic 
solutions.  The  dashed  black  lines  can  be  used  as  reference  points  for  comparing  different  figures.  Results  are  obtained  for  two  populations  with  size  A/1  =  500 
individuals  and  N2  =  1 000  individuals,  and  both  populations  are  initialized  with  a  single  infected  individual  l'0  =  Pa  =  1 . 
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for  the  stochastic  optimal  solution,  which  replicates  much  of  the  green  and  blue  high  probabili¬ 
ty  regions  above  the  intermediate  reference  line.  It  does  a  suboptimal  job  for  relatively  smaller 
epidemics  in  the  regions  where  the  amount  of  vaccine  ranges  from  V  =  400  to  V  =  1000. 

Fig.  10D  illustrates  the  corresponding  results  when  the  optimal  deterministic  protocol  is  ap¬ 
plied.  In  maximizing  P(E  <  Emax),  the  deterministic  protocol  underperforms  compared  to  the 
protocols  of  both  Fig.  10A  and  10C.  Comparatively,  the  deterministic  protocol  minimizes  the 
area  of  the  high  probability  (blue)  regions.  It  does  slightly  better  than  the  stochastic  optimum 
in  roughly  the  same  regions  where  the  stochastic  optimum  fails  compared  to  the  best  possible 
result,  from  about  V  =  400  to  V  =  750. 

This  shows  that  the  situation  does  indeed  become  more  complicated  when  one  looks  be¬ 
yond  optimizing  the  average  epidemic  size  ( E ) .  If  the  goal  is  to  keep  the  epidemic  below  some 
size,  given  some  amount  of  vaccine,  there  are  indeed  regions  where  the  deterministically 
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optimal  solution  will  yield  slightly  better  results.  Most  of  the  time  however,  optimizing  the  av¬ 
erage  stochastic  epidemic  size  gives  a  result  closer  to  the  best  possible  one  of  Fig.  10  A.  These 
figures  thus  indicate  that  the  average  epidemic  size  is  a  potentially  useful  metric  for  gauging  the 
effects  of  stochasticity  and  will  most  of  the  time  yield  a  solution  that  is  preferable  to  the 
deterministic  optimum. 


Discussion 

This  paper  illustrates  the  viability  and  power  of  developing  the  exact  numerical  solution  of  the 
master  equations,  done  here  for  the  stochastic  SIR  model.  We  modified  the  computational  al¬ 
gorithm  developed  by  Jenkinson  and  Goutsias  [15]  to  obtain  even  greater  numerical  efficiency. 
This  is  accomplished  by  removing  excess  states  of  the  system  which  have  no  probability  of 
occuring,  but  are  naturally  included  in  the  original  algorithm.  Even  more  significantly,  our 
work  and  that  of  Jenkinson  and  Goutsias  [15]  provide  proof  of  concept  for  obtaining  accurate, 
exact  solutions  for  SIR- type  models,  rather  than  relying  on  sampling  methods  [18]  or  approxi¬ 
mations  to  the  master  equation  [19]  [20].  Furthermore,  Black  and  Ross  recently  demonstrate 
an  efficient  and  numerically  stable  computational  methodology  that  computes  the  final  epi¬ 
demic  size  distribution,  for  a  broad  range  of  Markovian  SIR- type  models,  without  integration 
using  the  jump  chain  [21]. 

Our  analysis  focuses  on  the  fundamental  tradeoff  involving  ahocation  of  vaccine  between 
two  non-interacting  communities  of  different  size.  Our  procedure  involved  three  steps.  First, 
for  each  population  we  separately  calculate  the  probability  distribution  of  epidemic  sizes  for  a 
given  amount  of  vaccine.  Second,  we  evaluate  the  expected  epidemic  size  as  a  function  of  the 
amount  of  vaccine  in  each  population.  Third,  we  impose  a  constraint  on  the  total  amount  of 
vaccine  to  distribute  between  the  two  populations,  and  determine  the  optimal  allocation  which 
minimizes  the  expected  combined  epidemic  size  of  the  two  populations. 

We  obtain  several  results  that  serve  to  elaborate  and  refine  principles  first  identified  by  Kee¬ 
ling  and  Shattock  [2],  who  considered  the  corresponding  tradeoff  in  the  context  of  the  deter¬ 
ministic  SIR  model.  Where  the  deterministic  SIR  model  predicts  a  definite  epidemic  size  for 
any  given  set  of  parameters,  the  stochastic  SIR  model  produces  a  distribution,  the  characteris¬ 
tics  of  which  significantly  impact  protocols  for  optimal  allocation  of  vaccine.  Under  conditions 
that  promote  spread  of  the  epidemic  (i.e.,  the  reproductive  number  r0  >  1),  the  distribution  of 
epidemic  sizes  obtained  from  the  stochastic  SIR  model  is  bimodal  [14]  in  the  limit  of  large  pop¬ 
ulation  sizes,  consisting  of  a  peak  describing  terminal  infections,  that  fail  to  propagate  signifi¬ 
cantly  in  the  population,  and  a  peak  describing  large-scale  epidemics,  which  have  a  mean  size 
well  approximated  by  the  deterministic  size.  For  finite  population  sizes,  the  distinction  between 
terminal  infections  and  large-scale  epidemics  vanishes  at  a  value  of  r0  that  approaches  unity  as 
N  oo. 

Both  the  possibility  of  a  terminal  infection  and  the  width  of  the  distribution  of  the  large- 
scale  epidemic  sizes  contribute  significantly  to  differences  in  the  optimal  allocation  of  vaccine 
for  the  stochastic  model  compared  to  the  deterministic  case.  The  differences  are  most  signifi¬ 
cant  for  intermediate  ranges  of  vaccine.  In  contrast,  for  both  the  stochastic  and  deterministic 
cases,  when  vaccine  is  severely  limited  or  abundant,  there  is  little  or  no  difference  in  the  opti¬ 
mal  allocation  of  vaccine  between  the  two  models. 

Differences  in  optimal  allocations  are  amplified  for  intermediate  amounts  of  vaccine  be¬ 
cause  of  the  strong  switching  behavior  of  the  optimal  strategy.  This  switching  can  arise  in  both 
the  stochastic  and  deterministic  models,  but  at  different  points  quantitatively,  and  is  not  always 
observed  in  the  stochastic  case.  If  the  deterministic  protocol  is  applied  to  the  more  realistic  sto¬ 
chastic  description  of  the  epidemic  evolution  in  the  two  populations,  the  performance  is 
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suboptimal,  leading  to  a  greater  average  epidemic  size  than  would  occur  using  the  stochastic 
protocol.  The  difference  is  most  significant  for  smaller  values  of  r0  where  there  is  the  most  sig¬ 
nificant  probability  of  a  terminal  infection.  The  dependence  on  I0,  the  number  of  infected  indi¬ 
vidual,  is  more  complex  and  depends  on  r0,  but  in  the  limit  where  both  r0  and  I0  are  large,  the 
results  converge  to  those  of  the  deterministic  SIR  model.  In  the  absence  of  vaccine,  these  quan¬ 
tities  both  increase  the  relative  weight  in  the  peak  describing  terminal  infections. 

Keeling  and  Shattock  [2]  attribute  the  switching  behavior  to  the  property  of  herd  immunity, 
which  occurs  when  the  amount  of  vaccine  is  sufficient  to  prevent  the  epidemic  from  spreading 
significantly  in  the  population.  Herd  immunity  occurs  in  the  deterministic  SIR  model  when 
the  initial  effective  growth  rate  of  the  number  of  infected  individuals  in  the  population  becomes 
less  than  unity  [17].  While  the  optimal  deterministic  solution  approximately  distributes  vac¬ 
cine  in  a  manner  that  achieves  herd  immunity  in  the  largest  possible  population,  this  is  not  ex¬ 
actly  the  case.  More  precisely,  the  sharp  transitions  in  both  the  deterministic  and  stochastic 
models  arise  from  optimizing  the  overall  impact  of  the  vaccine  in  reducing  the  joint  epidemic 
size,  which  we  attribute  to  maximizing  the  overall  gain.  In  the  deterministic  model,  the  impact 
of  vaccine  on  epidemic  size  reduction  is  maximized  as  herd  immunity  is  approached.  For  the 
stochastic  model,  the  maximal  impact  typically  occurs  earlier,  and  in  some  cases,  there  is  no 
sharp,  intermediate  transition. 

Policies  involving  strong  switching  may  be  difficult  to  implement  publically,  as  one  commu¬ 
nity  could  be  reluctant  to  voluntarily  sacrifice  their  entire  vaccine  allocation  to  another  com¬ 
munity  in  favor  of  a  reduction  in  the  overall  epidemic  size.  In  contrast,  allocations  that  are 
proportional  to  population  size  are  likely  to  be  less  controversial  to  implement.  Our  results 
demonstrate  that  in  certain  scenarios  in  the  stochastic  model,  a  proportional  distribution  can 
be  justified  as  nearly  optimal.  For  intermediate  values  of  r0  (1  <  r0  SI  2.9)  and  small  initial  in¬ 
fected  populations,  we  find  that  switching  behavior  is  absent  in  the  optimal  stochastic  protocol. 
In  these  situations,  the  optimal  solution  is  reasonably  well  approximated  by  a  proportional  dis¬ 
tribution.  In  contrast,  in  situations  where  r0  and  I0  are  both  large,  and  large-scale  epidemics  are 
likely,  we  find  that  the  optimal  stochastic  protocol  is  more  similar  to  the  deterministic  case, 
and  proportional  distribution  results  in  significant  increases  in  the  overall  epidemic  size.  How¬ 
ever,  the  reduction  in  magnitude  of  the  gain  peaks  for  the  stochastic  model  in  Fig.  3,  compared 
to  the  deterministic  case,  indicate  that  the  overall  magnitude  of  the  benefit  (as  measured  by  re¬ 
duction  of  the  epidemic  size),  is  less  sensitive  to  the  precise  details  of  the  allocation  in  the  sto¬ 
chastic  model  than  it  is  in  the  corresponding  deterministic  case. 

Interestingly,  our  analysis  reveals  that  compared  to  the  deterministic  protocol,  the  stochastic 
protocol  that  minimizes  the  expected  epidemic  size,  also  overall  better  approximates  an  alter¬ 
native  target  based  on  specifying  a  maximum  tolerance  (or  threshold)  for  the  overall  epidemic 
size.  This  result  is  somewhat  surprising.  One  might  have  expected  the  deterministic  model  to 
be  more  accurate  in  this  case,  because  it  predicts  a  large-scale  epidemic  whenever  r0  >  1,  and  as 
such  might  have  captured  a  threshold  criterion  more  accurately.  The  fact  that  the  stochastic 
protocol  continues  to  outperform  the  deterministic  counterpart  provides  additional  impetus  to 
include  the  more  complete  and  accurate  stochastic  dynamics  of  epidemic  evolution  in 
further  studies. 

This  paper  isolates  the  tradeoff  in  vaccination  allocation  between  two  non-interacting  popu¬ 
lations,  prior  to  the  onset  of  widespread  disease,  in  order  to  illustrate  the  significance  of  the  full 
stochastic  solution  compared  to  deterministic  case.  Our  analysis  relies  on  some  strong  assump¬ 
tions,  particularly  the  assumption  of  non-interacting  populations.  The  extreme  switching  be¬ 
havior  in  the  deterministic  case  results  from  this  non-interaction.  It  is  less  clear  what  the 
optimal  policy  will  be  for  the  case  of  weakly  interacting  populations  in  the  stochastic  model. 
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One  might  speculate  that,  in  the  deterministic  limit,  the  presence  of  even  a  modest  amount  of 
interaction  yields  dynamics  of  a  single  population. 

Our  conceptual  framework  and  methods  can  potentially  be  generalized  to  include  increas¬ 
ingly  realistic  situations,  including  interacting  populations  and  real  time  allocation  of  vaccine 
as  the  epidemic  evolves.  In  these  scenarios,  we  anticipate  detailed  monitoring  of  stochastic  ef¬ 
fects,  as  well  as  incorporation  of  delays  associated  with  transportation  and  the  onset  of  immu¬ 
nity,  will  play  a  critical  role  in  determining  the  optimal  dynamic  protocol,  and  we  expect  that 
the  critical  differences  between  the  stochastic  and  deterministic  SIR  models  illustrated  here  will 
have  an  increasingly  significant  impact  in  identifying  protocols  that  aid  in  minimizing  the 
overall  epidemic  size. 

The  hope  is  that  the  systematic  study  of  such  tradeoffs  will  shed  light  on  the  development  of 
effective  policies.  For  example,  in  the  case  of  epidemic  outbreak  in  a  localized  geographic  re¬ 
gion,  government  officials  might  have  to  decide  whether  to  allocate  scare  vaccination  doses  ex¬ 
clusively  to  that  region  or  to  allocate  the  vaccine  proportionately  for  population  as  a  whole.  In 
situations  where  vaccine  doses  have  been  prepositioned  geographically,  the  question  of  “giving 
away”  vaccines  from  one  region  to  another  will  be  of  intense  debate.  Thus,  issues  of  fairness 
will  complicate  decisions  even  more.  Identifying  policies  that  are  close-to-optimal  and  can  ac¬ 
tually  be  implemented  is  an  important  topic  for  future  research. 

Acknowledgments 

The  authors  thank  Charles  Lieou  and  Kimberly  Schlesinger  for  helpful  discussions  and  feed¬ 
back.  This  work  was  supported  by  an  Office  of  Naval  Research  MURI  Grant  No. 

DMR0606092,  the  David  and  Lucile  Packard  Foundation,  the  Stansberry  Fellowship  through 
the  CCS  SURF  foundation,  and  the  Institute  for  Collaborative  Biotechnologies  through  con¬ 
tract  no.  W911NF-09-D-0001  from  the  U.S.  Army  Research  Office.  The  funders  had  no  role  in 
study  design,  data  collection  and  analysis,  decision  to  publish,  or  preparation  of 
the  manuscript. 

Author  Contributions 

Conceived  and  designed  the  experiments:  EY  SS  DA  JC.  Performed  the  experiments:  EY  SS  JC. 
Analyzed  the  data:  EY  JC.  Contributed  reagents/materials/analysis  tools:  EY  SS  JC.  Wrote  the 
paper:  EY  DA  JC. 

References 

1.  Leroux-Roels  I,  Leroux-Roels  G  (2009)  Current  status  and  progress  of  prepandemic  and  pandemic  in¬ 
fluenza  vaccine  development.  Expert  Review  of  Vaccines  8:  401  -423. 

2.  Keeling  MJ,  Shattock  A  (2012)  Optimal  but  unequitable  prophylactic  distribution  of  vaccine.  Epidemics 
4:  78-85.  doi:  10.1016/j.epidem.2012.03.001  PMID:  22664066 

3.  Wu  JT,  Riley  S,  Leung  GM  (2007)  Spatial  considerations  for  the  allocation  of  pre-pandemic  inuenza 
vaccination  in  the  united  states.  Proc  R  Soc  B  274:  281 1-2817.  doi:  10.1098/rspb.2007.0893  PMID: 
17785273 

4.  Knipl  DH,  Rost  G  (201 1 )  Modelling  the  strategies  for  age  specific  vaccination  scheduling  during  influen¬ 
za  pandemic  outbreaks.  Mathematical  Biosciences  and  Engineering  8: 123-139.  doi:  10.3934/mbe. 
2011.8. 123  PMID:  21361404 

5.  Kermack  WO,  McKendrick  AG  (1 927)  A  contribution  to  the  mathematical  theory  of  epidemics.  Proc  R 
SocLondA  1 15:  700-721.  doi:  10.1098/rspa.  1927.01 18 

6.  Anderson  M,  May  R  (1 982)  Directly  transmitted  infections  diseases:  Control  by  vaccination.  Science 
215: 10531060.  doi:  10.1 126/science.7063839 

7.  Gaff  H,  Schaefer  E  (2009)  Optimal  control  applied  to  vaccination  and  treatment  strategies  for  various 
epidemiological  models.  Mathematical  Biosciences  and  Engineering  6:  469-492.  doi:  10.3934/mbe. 
2009.6.469  PMID:  19566121 


PLOS  ONE  |  DOI:  10. 1371 /journal. pone. 01 1 5826  February  1 7,  201 5 


24/25 


PLOS 


ONE 


Optimal  Vaccination  in  a  Stochastic  Epidemic  Model 


8.  Klepac  P,  Laxminarayan  R,  Grenfell  BT  (2011)  Synthesizing  epidemiological  and  economic  optima  for 
control  of  immunizing  infections.  Proc  Natl  Acad  Sci  U  S  A  108: 14366-14370.  doi:  10.1073/pnas. 

1 1 01 6941 08  PMID:  218251 29 

9.  Tildesley  MJ,  Savill  NJ,  Shaw  DJ,  Deardon  R,  Brooks  SP,  et  al.  (2006)  Optimal  reactive  vaccination 
strategies  for  an  outbreak  of  Foot-and-mouth  disease  in  Great  Britain.  Nature  440:  83-86.  doi:  1 0. 

1038/nature04324  PMID:  16511494 

10.  Keeling  MJ,  Woolhouse  MEJ,  Shaw  DJ,  Matthews  L,  Chase-Topping  M,  et  al.  (2001)  Dynamics  of  the 
2001  UK  Foot-and-mouth  epidemic:  Stochastic  dispersal  in  a  heterogeneous  landscape.  Science  294: 
813-817.  doi:  10.1 126/science.  1065973  PMID:  11679661 

1 1 .  Groenendaal  H,  Nielen  M,  Jalvingh  AW,  Horst  SH,  Galligan  DT,  et  al.  (2002)  A  simulation  of  Johne’s 
disease  controls.  Prev  Vet  Med  54:  225-245.  doi:  10. 101 6/S01 67-5877(02)00027-2  PMID:  12114011 

12.  Grenfell  BT,  Bjornstad  ON,  Kappey  J  (2001)  Travelling  waves  and  spatial  hierarchies  in  measles  epi¬ 
demics.  Science  414:  716-723. 

13.  Hufnagel  L,  Brockmann  D,  GeiselT  (2003)  Forecast  and  control  of  epidemics  in  a  globalized  world. 
Proc  Natl  Acad  Sci  U  S  A  101 : 15124-15129.  doi:  10.1073/pnas.0308344101 

1 4.  Gordillo  LF,  Marion  SA,  Martin-Lof  A,  Greenwood  PE  (2008)  Bimodal  epidemic  size  distributions  for 
near-critical  sir  with  vaccination.  Bulletin  of  Mathematical  Biology  70:  589-602.  doi:  10.1007/si  1538- 
007-9269-y  PMID:  17992563 

15.  Jenkinson  G,  Goutsias  J  (2012)  Numerical  integration  of  the  master  equation  in  some  models  of  sto¬ 
chastic  epidemiology.  PLoSONE  7:  e36160.  doi:  10. 1371/journal. pone. 0036160  PMID:  22615755 

1 6.  Lalley  SP,  Perkins  EA,  Zheng  X  (2014)  A  phase  transition  for  measure-valued  sir  epidemic  processes. 
The  Annals  of  Probability  42:  237-310.  doi:  10.1214/13-AOP846 

17.  De  Jong  MC,  Bouma  A  (2001)  Herd  immunity  after  vaccination:  Howto  quantify  it  and  howto  use  it  to 
halt  disease.  Vaccine  19:  2722-2728.  doi:  10.1016/S0264-410X(00)00509-0  PMID:  1 1257415 

18.  Petrovic  N,  Alderson  DL,  Carlson  JM  (2012)  Dynamic  resource  allocation  in  disaster  response:  Trade¬ 
offs  in  wildfire  suppression.  PLoS  ONE  7:  e33285.  doi:  10. 1371/journal. pone. 0033285  PMID: 
22514605 

1 9.  Black  A,  McKane  A  (201 0)  Stochasticity  in  staged  models  of  epidemics:  quantifying  the  dynamics  of 
whooping  cough.  J  RSoc  Interface  7:  1219-1227.  doi:  10. 1098/rsif. 2009.051 4  PMID:  20164086 

20.  Munsky  B,  Khammash  M  (2006)  The  finite  projection  algorithm  for  the  solution  of  the  chemical  master 
equation.  J  Chem  Phys  124:  044104.  doi:  10.1063/1 .2145882  PMID:  16460146 

21 .  Black  AJ,  Ross  JV  (201 4)  Computation  of  epidemic  final  size  distributions.  ArXiv  e-prints. 


PLOS  ONE  |  DOI:  10. 1371 /journal. pone. 01 1 5826  February  1 7,  201 5 


25/25 


